#Libraries
##### Begin Block 1
## John fix
setwd("~/AnimeshReview/LauraPaper/Journal of Translational Medicine /LauraFinal/")
#load("finalcommentsplotsandall.RData")
library(readxl)
library(data.table)
library(tidyverse)
library(dplyr)
library(devtools)
library(ggcorrplot)
library(car)
library(ggpubr)
library(glmnet)
library(summarytools)
library(knitr)
library(htmltools)
library(corrplot)
library(caret)
library(factoextra)
library(Metrics)
library(readr)
library(gplots)
library(dplyr)
library(stringr)
library(readxl)
library(plotly)
library(e1071)
library(ggplot2)
library(reshape2)
library(multtest)
library(ROCR)
library(gridExtra)
library(MLmetrics)
#### Begin Load Laura's Functions
doubleAUCfun <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
model3<-glm(Outcome~Predictor+Predictor2,data=yy3, family=binomial(link='logit'))
new3<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
}
doubleAUCfunNB <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
model4<-naiveBayes(Outcome~Predictor+Predictor2,data=yy3)
new4<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p4<-predict(model4,new4,type="raw")
#without raw write confusion matrix directly
prob2<-NULL
for (i in 1:dim(p4)[1]){
prob2[i]<-p4[i,2]/p4[i,1]
}
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(prob2, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
}
doubleAUCfunRFCross <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
control <- trainControl(method="repeatedcv", number=10, repeats=3) #Is this the only thing needed for cv
seed <- 7
metric <- "Accuracy"
set.seed(seed)
#mtry <- 2
tunegrid <- expand.grid(.mtry=c(1:6))
set.seed(seed)
rf_default <- train(Outcome~Predictor+Predictor2,data=yy3, method="rf", metric=metric, tuneGrid=tunegrid, trControl=control)
#print(rf_default)
new4<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
pred=predict(rf_default, new4,type="prob")
pred2=predict(rf_default, new4)
prob2<-NULL
for (i in 1:dim(pred)[1]){
prob2[i]<-pred[i,2]/pred[i,1]
}
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(prob2, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
}
doubleAUCfunSVM <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
seed <- 7
set.seed(seed)
model3<-svm(Outcome~Predictor+Predictor2,data=yy3)
new3<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p3<-predict(model3,new3,decision.values = TRUE)
p4<-attr(p3,"decision.values")
new5<-data.frame(Outcome=ytest)
pr3<- prediction(p4, new5)
#table(p3, new5)
#confusionMatrix(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
}
doubleAUCfunSVMCross <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
seed <- 7
set.seed(seed)
model4<-svm(Outcome~Predictor+Predictor2,data=yy3,method="C-classification",
kernel="radial", gamma = 0.01, cost = 100,cross=10, probability=TRUE)
new4<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p4<-predict(model4,new4,decision.values = TRUE)
p5<-attr(p4,"decision.values")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p5, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
}
doublePlusfun <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
model3<-glm(Outcome~Predictor+Predictor2,data=yy3, family=binomial(link='logit'))
new3<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr <- prediction(p3, new5)
Acc <- performance(pr, measure="acc")
AccV<-Acc@y.values[[1]][max(which(Acc@x.values[[1]] >= 0.5))]
Sens <- performance(pr, measure= "sens")
SensV<-Sens@y.values[[1]][max(which(Sens@x.values[[1]] >= 0.5))]
Spec <- performance(pr, measure= "spec")
SpecV<-Spec@y.values[[1]][max(which(Spec@x.values[[1]] >= 0.5))]
Prec <- performance(pr, measure= "prec")
PrecV<-Prec@y.values[[1]][max(which(Prec@x.values[[1]] >= 0.5))]
AllV<-data.frame(Vector=c(AccV, SensV, SpecV, PrecV))
return(AllV)
}
doubleROCfun <- function(xtrain, ytrain,xtest,ytest,s,k) {
yy3<-data.frame(Outcome=ytrain, Predictor=xtrain[,s],Predictor2=xtrain[,k])
model3<-glm(Outcome~Predictor+Predictor2,data=yy3, family=binomial(link='logit'))
new3<-data.frame(Predictor=xtest[,s],Predictor2=xtest[,k])
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
result<-prf3
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
return(result)
}
MeansNames <- function(doubleAUCR,trial){
doubleAUCR<-as.data.frame(doubleAUCR)
doubleAUCR<-mutate(doubleAUCR, Means=rowMeans(doubleAUCR))
row.names(doubleAUCR)<-trial
return(doubleAUCR)
}
multipleAUCfun <- function(xtrain,ytrain,xtest,ytest) {
yy3<-data.frame(Outcome=ytrain,xtrain)
model3<-glm(Outcome~.,data=yy3, family=binomial(link='logit'))
new3<-data.frame(xtest)
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
return(result2)
}
multipleAUCfunNB <- function(xtrain, ytrain,xtest,ytest) {
yy3<-data.frame(Outcome=ytrain,xtrain)
model4<-naiveBayes(Outcome~.,data=yy3)
new3<-data.frame(xtest)
p4<-predict(model4,new3,type="raw")
#without raw write confusion matrix directly
prob2<-NULL
for (i in 1:dim(p4)[1]){
prob2[i]<-p4[i,2]/p4[i,1]
}
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(prob2, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
return(result2)
}
multipleROCfun <- function(xtrain,ytrain,xtest,ytest) {
yy3<-data.frame(Outcome=ytrain,xtrain)
model3<-glm(Outcome~.,data=yy3, family=binomial(link='logit'))
new3<-data.frame(xtest)
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
return(prf3)
}
multipleROCfunNB <- function(xtrain, ytrain,xtest,ytest) {
yy3<-data.frame(Outcome=ytrain,xtrain)
model4<-naiveBayes(Outcome~.,data=yy3)
new3<-data.frame(xtest)
p4<-predict(model4,new3,type="raw")
#without raw write confusion matrix directly
prob2<-NULL
for (i in 1:dim(p4)[1]){
prob2[i]<-p4[i,2]/p4[i,1]
}
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(prob2, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
return(prf3)
}
singleAUCfun <- function(xtrain, ytrain,xtest,ytest,s) {
yy2<-data.frame(Outcome=ytrain, Predictor=xtrain[,s])
model2<-glm(Outcome~Predictor,data=yy2,family=binomial(link='logit'))
new2<-data.frame(Predictor=xtest[,s])
p<-predict(model2,new2,type = "response")
new4<-data.frame(Outcome=ytest)
pr <- prediction(p, new4)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
auc <- performance(pr, measure = "auc")
auc <- auc@y.values[[1]]
return(auc)
}
singlePlusfun <- function(xtrain, ytrain,xtest,ytest,s) {
yy2<-data.frame(Outcome=ytrain, Predictor=xtrain[,s])
model2<-glm(Outcome~Predictor,data=yy2,family=binomial(link='logit'))
new2<-data.frame(Predictor=xtest[,s])
p<-predict(model2,new2,type = "response")
new4<-data.frame(Outcome=ytest)
pr <- prediction(p, new4)
Acc <- performance(pr, measure="acc")
AccV<-Acc@y.values[[1]][max(which(Acc@x.values[[1]] >= 0.5))]
Sens <- performance(pr, measure= "sens")
SensV<-Sens@y.values[[1]][max(which(Sens@x.values[[1]] >= 0.5))]
Spec <- performance(pr, measure= "spec")
SpecV<-Spec@y.values[[1]][max(which(Spec@x.values[[1]] >= 0.5))]
Prec <- performance(pr, measure= "prec")
PrecV<-Prec@y.values[[1]][max(which(Prec@x.values[[1]] >= 0.5))]
AllV<-data.frame(Vector=c(AccV, SensV, SpecV, PrecV))
return(AllV)
}
singleROCfun <- function(xtrain, ytrain,xtest,ytest,s) {
yy2<-data.frame(Outcome=ytrain, Predictor=xtrain[,s])
model2<-glm(Outcome~Predictor,data=yy2,family=binomial(link='logit'))
new2<-data.frame(Predictor=xtest[,s])
p<-predict(model2,new2,type = "response")
new4<-data.frame(Outcome=ytest)
pr <- prediction(p, new4)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
resultU2<-prf
tomeany<-as.data.frame(prf@y.values)
tomeanx<-as.data.frame(prf@x.values)
auc <- performance(pr, measure = "auc")
auc <- auc@y.values[[1]]
return(resultU2)
}
#### End Load Laura's Functions, edit with seeds for SVM and RF cuz random
set.seed(132)
#LASSO ####
#setwd("../R code /")
# Hazeldine_1h <- read_csv("CF_ISS1.csv")
# savee<-as.factor(Hazeldine_1h$Label)
# names(Hazeldine_1h)<-sapply(1:dim(Hazeldine_1h)[2], function(i){paste0("First ",names(Hazeldine_1h)[i])})
# Hazeldine_1h<-Hazeldine_1h[,-(1:3)]
# Hazeldine_1h<-Hazeldine_1h[,-(77:78)]
#
# Hazeldine_2h <- read_csv("CF_ISS2.csv")
# names(Hazeldine_2h)<-sapply(1:dim(Hazeldine_2h)[2], function(i){paste0("Second ",names(Hazeldine_2h)[i])})
# Hazeldine_2h<-Hazeldine_2h[,-(1:3)]
# Hazeldine_2h<-Hazeldine_2h[,-(79:80)]
#
# Hazeldine_3h <- read_csv("CF_ISS3.csv")
# Hazeldine_3h<-Hazeldine_3h[,-(1:3)]
# names(Hazeldine_3h)<-sapply(1:dim(Hazeldine_3h)[2], function(i){paste0("Third ",names(Hazeldine_3h)[i])})
#
#
# #The features from all three data points are merged together, as same name, add descriptor to identify each time point.
# #Label accounts for the MODS development
# All2<-data.frame(Hazeldine_1h,Hazeldine_2h,Hazeldine_3h,Label=savee)
#
#
# saveRDS(All2,"All2.Laura.rds")
# saveRDS(savee,"savee.Laura.rds")
### Load Laura's Data
#setwd("/Users/j.williams/AnimeshReview/LauraPaper/Journal of Translational Medicine /LauraFinal/Sept18_LauraRes")
#All2 <- readRDS("All2.Laura.rds"); savee <- readRDS("savee.Laura.rds")
### Load AML Data
#setwd("/Users/j.williams/AnimeshReview/LauraPaper/Journal of Translational Medicine /LauraFinal/Sept18_AMLRes")
#All2 <- readRDS("All2.AML.rds"); savee <- readRDS("savee.AML.rds")
### Load Pancreatic Cancer Dta
setwd("/Users/j.williams/AnimeshReview/LauraPaper/Journal of Translational Medicine /LauraFinal/Sept18_PancreaticRes")
All2 <- readRDS("All2.Pancreatic.rds"); savee <- readRDS("savee.Pancreatic.rds")
#TRAIN
### set n to 10 for now to make quicker
n<-100 #Number of LASSO and EN loops. Per loop a model with features selected. The 80% most popular features (those that appear in 40/50 or more) will be selected
#N and nn control the number of models created for combinatorial analysis to help in further feature selection - Highest AUC values chosen.
N <-10
nn <-25
ErrorsFinEN<-vector(mode="double", length=n)
BetasFinEN<-vector(mode="character", length=n)
LambdaFinEN<-vector(mode="double", length=n)
BNumFinEN<-vector(mode="double", length=n)
see2EN<-data.frame(All="All")
LauCoef1<-data.frame(Coeff="See",stringsAsFactors=FALSE)
BetasTodo<-data.frame(Features="Name",Coefficients=1)
ListError<-vector(mode="double", length=n)
BetasFin<-vector(mode="character", length=n)
LambdaFin<-vector(mode="double", length=n)
BNumFin<-vector(mode="double", length=n)
see2<-data.frame(All="All")
LauCoef1L<-data.frame(Coeff="See",stringsAsFactors=FALSE)
BetasTodoL<-data.frame(Features="Name",Coefficients=1)
##### End Block 1
##### Begin Block 2
for (i in 1:n){
smp_size = floor(0.75 * nrow(All2))
#train_ind = sample(seq_len(nrow(All2)), size = smp_size)
train_ind <- caret::createDataPartition(All2$Label, p = 0.75, list = FALSE, times = 1)
#Training set
train = All2[train_ind, ]
#Test set
test = All2[-train_ind, ]
#Creates matrices for independent and dependent variables.
## try something like
####### xtrain <- train[ , !(names(train) %in% "Label")]
xtrain <- train[ , !(names(train) %in% "Label")] %>% as.matrix()
#xtrain = model.matrix(Label~. -1, data = train)
ytrain = train$Label
xtest <- test[ , !(names(test) %in% "Label")] %>% as.matrix()
#xtest = model.matrix(Label~. -1, data = test)
ytest = test$Label
#Choose lambda value that minimize missclassification error.
#0.5 as elastic nets, all variables with EN are based on ElasticNets analysis. 100 lambdas sampled with 10 cross validation for each, already internalized in method
CVEN=cv.glmnet(xtrain,ytrain,family="binomial",type.measure="class",alpha=0.5,nlambda=100)
attach(CVEN)
Lambda.BestEN<-CVEN$lambda.min #can be either minimum or 1 standard deviation
print(Lambda.BestEN)
CVFinEN=glmnet(xtrain,ytrain,family="binomial",alpha=0.5,lambda=Lambda.BestEN)
CoefEN<-coef(CVFinEN) #Beta coefficients obtained from here
InterceptEN<-CoefEN@x[1]
BetasEN<-CVFinEN$beta
Betas2EN<-data.frame(Features=BetasEN@Dimnames[[1]][BetasEN@i+1], Coefficients=BetasEN@x) #Beta coefficients names stored here
CVPred1EN = predict(CVFinEN, family="binomial", s=Lambda.BestEN, newx = xtest,type="class") #predict in test set to obtain confusion matrix
#Calculate error for categorical values
ytest2<-as.factor(ytest)
ResultsEN<-table(CVPred1EN,ytest)
confusionMatrix(as.factor(CVPred1EN),ytest)
AccuracyEN<-(ResultsEN[1]+ResultsEN[4])/sum(ResultsEN[1:4])
ErrorEN<-1-AccuracyEN
LauCoef<-Betas2EN$Coefficients
LauCoefEN<-data.frame(Coeff=LauCoef,stringsAsFactors=FALSE)
LauCoef1<-rbind(LauCoef1,LauCoefEN)
BetasTodo<-rbind(BetasTodo,Betas2EN) #store coefficients and store betas
seeEN<-Betas2EN$Features
seeEN1<-data.frame(All=seeEN)
see2EN<-rbind(see2EN,seeEN1) #all beta names stored
mEN<-count(see2EN, All) #frequency of the betas stored counted
see3EN<-toString(seeEN)
ErrorsFinEN[i]<-ErrorEN #error of the model stored
BetasFinEN[i]<-see3EN #name of features the model used
BNumFinEN[i]<-length(seeEN) #number of features the model used
LambdaFinEN[i]<-Lambda.BestEN #lambda chosen for model
detach(CVEN)
#Change between Lasso and EN, alpha=1 (*)
CV=cv.glmnet(xtrain,ytrain,family="binomial",type.measure="class",alpha=1,nlambda=100)
attach(CV)
Lambda.Best<-CV$lambda.min
CVFin=glmnet(xtrain,ytrain,family="binomial",alpha=1,lambda=Lambda.Best)
Coef<-coef(CVFin)
Intercept<-Coef@x[1]
Betas<-CVFin$beta
Betas2<-data.frame(Features=Betas@Dimnames[[1]][Betas@i+1], Coefficients=Betas@x)
CVPred1 = predict(CVFin, family="binomial", s=Lambda.Best, newx = xtest,type="class")
#Calculate error for categorical values
ytest2<-as.factor(ytest)
confusionMatrix(as.factor(CVPred1),ytest)
Results<-table(CVPred1,ytest)
Accuracy<-(Results[1]+Results[4])/sum(Results[1:4])
Error<-1-Accuracy
#visual display of for
BetasTodoL<-rbind(BetasTodoL,Betas2)
see<-Betas2$Features
see1<-data.frame(All=see)
see2<-rbind(see2,see1)
m<-count(see2, All)
see3<-toString(see)
ListError[i]<-Error
BetasFin[i]<-see3
BNumFin[i]<-length(see)
LambdaFin[i]<-Lambda.Best
detach(CV)
}
## [1] 0.6689277
## [1] 0.6303686
## [1] 0.7294298
## [1] 0.2883708
## [1] 0.4733006
## [1] 0.5943159
## [1] 0.235993
## [1] 0.6584265
## [1] 0.1317309
## [1] 0.395802
## [1] 0.6748829
## [1] 0.6907273
## [1] 0.4808845
## [1] 0.2339522
## [1] 0.6926208
## [1] 0.4017377
## [1] 0.6815938
## [1] 0.3682275
## [1] 0.6878564
## [1] 0.5334466
## [1] 0.25522
## [1] 0.7062928
## [1] 0.5021274
## [1] 0.2454332
## [1] 0.3415769
## [1] 0.1445318
## [1] 0.2492792
## [1] 0.7286193
## [1] 0.4040431
## [1] 0.5608647
## [1] 0.4575239
## [1] 0.639378
## [1] 0.1181126
## [1] 0.194737
## [1] 0.1059382
## [1] 0.4479945
## [1] 0.3634272
## [1] 0.4787447
## [1] 0.5555517
## [1] 0.5364108
## [1] 0.1343051
## [1] 0.4800724
## [1] 0.3865033
## [1] 0.6247686
## [1] 0.2318932
## [1] 0.1884252
## [1] 0.6696531
## [1] 0.5132077
## [1] 0.502099
## [1] 0.629618
## [1] 0.06499661
## [1] 0.6696552
## [1] 0.3857286
## [1] 0.2742923
## [1] 0.2977953
## [1] 0.607225
## [1] 0.2129634
## [1] 0.4963329
## [1] 0.4782131
## [1] 0.2843202
## [1] 0.3724457
## [1] 0.4973665
## [1] 0.6293572
## [1] 0.251071
## [1] 0.7596586
## [1] 0.4509715
## [1] 0.4194573
## [1] 0.2410539
## [1] 0.6456791
## [1] 0.7355681
## [1] 0.631905
## [1] 0.69536
## [1] 0.5778142
## [1] 0.4760281
## [1] 0.642023
## [1] 0.36937
## [1] 0.4357793
## [1] 0.2766556
## [1] 0.5132012
## [1] 0.1636617
## [1] 0.7393992
## [1] 0.4958403
## [1] 0.172681
## [1] 0.6463179
## [1] 0.5903157
## [1] 0.393137
## [1] 0.7239418
## [1] 0.3621122
## [1] 0.6498011
## [1] 0.2339883
## [1] 0.5483706
## [1] 0.2081683
## [1] 0.5177668
## [1] 0.6192567
## [1] 0.609796
## [1] 0.2968183
## [1] 0.6362316
## [1] 0.4532429
## [1] 0.6682868
## [1] 0.4565895
#Visualizing data from LASSO and EN ####
#obtain in a data frame all error, betas names, number and lamda for the N models for each lasso and EN
All_info<-data.frame(Error=ListError, BetasNames=BetasFin, BetasNum=BNumFin, Lambda=LambdaFin)
All_infoEN<-data.frame(Error=ErrorsFinEN, BetasNames=BetasFinEN, BetasNum=BNumFinEN, Lambda=LambdaFinEN)
m<-m[-1,]
mEN<-mEN[-1,]
Final_LASSO<-m[order(-m$n),] #order highest frequencies above and filter with those that appear more than 80% of times
Final_LASSO1<-filter(Final_LASSO,n>40) #threshold selected - 80%
Final_EN<-mEN[order(-mEN$n),]
Final_EN1<-filter(Final_EN,n>40)
Final_Plot_Names<-filter(Final_EN,n>40)
outputVenn2<-venn(list(EN= Final_EN$All, LASSO = Final_LASSO$All))
outputVenn<-venn(list(EN= Final_EN1$All, LASSO = Final_LASSO1$All))
Freqs<-m[order(-m$n),]
num<-length(Freqs$All)
Freqs$All <- factor(Freqs$All, levels = Freqs$All[order(-Freqs$n)]) #plot in a bar graph the frequencies of ocurrance of the betas and order from highest to smalles
ggplot(Freqs, aes(All, n))+geom_bar(stat="identity")+theme(axis.text.x = element_text(size=8, angle=90))+ggtitle("LASSO features")
FreqsEN<-mEN[order(-mEN$n),]
numEN<-length(FreqsEN$All)
FreqsEN$All <- factor(FreqsEN$All, levels = FreqsEN$All[order(-FreqsEN$n)])
ggplot(FreqsEN, aes(All, n))+geom_bar(stat="identity")+theme(axis.text.x = element_text(size=8, angle=90))+ggtitle("EN features")
#plot of how many times each feature appears. Most important will appear in all models = N.
#Boxplot with Betas and its coefficients
Boxplot1<-BetasTodo[BetasTodo$Features %in% Final_EN1$All,] #see which features appear in the filtered features (80%) and obtain the coefficients associated.
ggplot(Boxplot1,aes(Boxplot1$Features,Boxplot1$Coefficients))+geom_boxplot()+geom_jitter()
Boxplot1["Method"]<-as.factor("EN")
Boxplot2<-BetasTodoL[BetasTodoL$Features %in% Final_LASSO1$All,]
ggplot(Boxplot2,aes(Boxplot2$Features,Boxplot2$Coefficients))+geom_boxplot()+geom_jitter()
Boxplot2["Method"]<-as.factor("LASSO")
Fin_Boxplot<-rbind(Boxplot1,Boxplot2) #Unite both boxplots LASSO and EN
ggplot(Fin_Boxplot,aes(Fin_Boxplot$Features,Fin_Boxplot$Coefficients))+geom_boxplot(aes(color=Method))+geom_jitter()+theme(axis.text.x = element_text(size=12, angle=90))+ggtitle("Beta coefficients EN and LASSO")+ labs(x = "Features",y="Beta Coefficients")
All_Feat<-rbind(Final_LASSO1,Final_EN1)
All_Feat2<-unique(All_Feat$All) #select the filtered betas in both.
#All_Feat2<-data.frame(All_Feat2)
Betas_select<-All2[,colnames(All2[,intersect(gsub("`", "", All_Feat2), colnames(All2))])]
Betas_select["Label"]<-All2$Label
Betas_select<-data.frame(Betas_select)
#Beta_select are the final features selected.
##### End Block 2
##### Begin Block 3
###Combinatorial analysis ####
NumVar<-length(Betas_select)
maxCombF<-data.frame(Name=toString("Trial"),AUC=0)
auc3max<-0
cont<-0
cont2<-0
for (p in 1:N){ #N different measurements of AUC values, mean done at the end.
smp_size<-floor(0.65 * nrow(Betas_select))
#set.seed(907)
# train_ind<-sample(seq_len(nrow(Betas_select)), size = smp_size)
train_ind <- caret::createDataPartition(All2$Label, p = 0.65, list = FALSE, times = 1)
# Training set
train<-Betas_select[train_ind, ]
# Test set
test<-Betas_select[-train_ind, ]
xtrain <- train[ , !(names(train) %in% "Label")] %>% as.matrix()
xtest <- test[ , !(names(test) %in% "Label")] %>% as.matrix()
# xtrain<-model.matrix(Label~. -1, data = train)
ytrain<-train$Label
# xtest<-model.matrix(Label~. -1, data = test)
ytest<-test$Label
xtest<-data.frame(xtest)
xtrain<-data.frame(xtrain)
y<-Betas_select[,NumVar]
X<-Betas_select[,1:(NumVar-1)]
levels(ytrain)[1]<-"0"
levels(ytrain)[2]<-"1"
levels(ytest)[1]<-"0"
levels(ytest)[2]<-"1"
for (k in 1:nn){
columns<-c(1:dim(xtrain)[2])
columns<-sample(columns)
d<-xtrain[,columns]
for (i in 1:dim(xtrain)[2]){
for (j in 1:dim(xtrain)[2]){
yy3<-data.frame(Outcome=ytrain,d[i:j])
model3<-glm(Outcome~.,data=yy3, family=binomial(link='logit'))
dd<-xtest[,columns]
new3<-data.frame(dd[i:j])
p3<-predict(model3,new3,type="response")
new5<-data.frame(Outcome=ytest)
pr3 <- prediction(p3, new5)
prf3 <- performance(pr3, measure = "tpr", x.measure = "fpr")
auc3 <- performance(pr3, measure = "auc")
auc3 <- auc3@y.values[[1]]
result2<-auc3
if (auc3>auc3max){
if (i>j){
maxComb<-data.frame(Name=toString(names(d)[j:i]),AUC=auc3)
auc3max<-auc3
cont<-cont+1
}
else{
maxComb<-data.frame(Name=toString(names(d)[i:j]),AUC=auc3)
auc3max<-auc3
cont<-cont+1
}
}
else{
cont2<-cont2+1
}
}
}
}
maxCombF<-rbind(maxCombF,maxComb)
auc3max<-0
}
maxCombF2<-maxCombF[order(maxCombF$AUC),]
names<-maxCombF2$Name[dim(maxCombF2)[1]]
names1<-as.character(names)
names1<-strsplit(names1,", ")
names<-as.data.frame(names1)
Betas_select2<-All2[,colnames(All2[,intersect(gsub("`", "", names[,1]), colnames(All2))])]
Betas_select2["Label"]<-savee
Betas_select2<-as.data.frame(Betas_select2)
##### End Block 3
##### Begin Block 4
#Random mix between variables and AUC value obtained through GLM model, rough approximation of future performance. Could substitute/add to filtering with clinical knowledge.
#Betas_select2, final features selected.
NumVar<-length(Betas_select2)
#ISSs<-c("Third.NISS","Third.ISS")
#We would obtain ISS or NISS as important features, never together as they are highly correlated. Have to know in which position of the new selected features it is found given we perfome a bivariate model with it and added feature.
#where<-match(colnames(Betas_select2),ISSs)
#for (u in 1:dim(Betas_select2)[2]){
# if (!(is.na(where[u]))){
# if (where[u]==1){
# CompareISS<-u
# NISShere<-1
# }else{
# CompareISS<-u
# NISShere<-0
# }
# }
#}
N<-1000 #number of models produced (both permuted (random) and real) - for all univariate (Each feature), bivariate and multivariate.
#Real program runs with 1000 - Takes 12 hours.
multipleAUC<-matrix(rnorm(2),1,N)
multipleAUCR<-matrix(rnorm(2),1,N)
multipleAUCNB<-matrix(rnorm(2),1,N)
multipleAUCNBR<-matrix(rnorm(2),1,N)
multipleROC<-matrix(as.list(rnorm(2)),1,N)
multipleROCR<-matrix(as.list(rnorm(2)),1,N)
multipleNBROC<-matrix(as.list(rnorm(2)),1,N)
multipleNBROCR<-matrix(as.list(rnorm(2)),1,N)
singleROC<-list()
doubleROC<-list()
singleROCR<-list()
doubleROCR<-list()
doublePlus<-list()
singlePlus<-list()
singleAUC<-matrix(rnorm(2),NumVar-1,N)
doubleAUC<-matrix(rnorm(2),(NumVar-1),N)
singleAUCR<-matrix(rnorm(2),NumVar-1,N)
doubleAUCR<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCSVMR<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCSVM<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCSVMCrossR<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCSVMCross<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCRFCross<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCRFCrossR<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCNBR<-matrix(rnorm(2),(NumVar-1),N)
doubleAUCNB<-matrix(rnorm(2),(NumVar-1),N)
MatsingleROC<-matrix(as.list(rnorm(2)),NumVar-1,N)
MatsingleROCR<-matrix(as.list(rnorm(2)),NumVar-1,N)
MatdoubleROC<-matrix(as.list(rnorm(2)),(NumVar-1),N)
MatdoubleROCR<-matrix(as.list(rnorm(2)),(NumVar-1),N)
MatsinglePlus<-matrix(as.list(rnorm(2)),NumVar-1,N)
MatdoublePlus<-matrix(as.list(rnorm(2)),(NumVar-1),N)
MatsinglePlusR<-matrix(as.list(rnorm(2)),NumVar-1,N)
MatdoublePlusR<-matrix(as.list(rnorm(2)),(NumVar-1),N)
##### End Block 4
##### Begin Block 5
for (j in 1:N){ #N different measurements of AUC values, mean done at the end.
smp_size<-floor(0.65 * nrow(Betas_select2))
#set.seed(907)
#train_ind<-sample(seq_len(nrow(Betas_select2)), size = smp_size)
train_ind <- caret::createDataPartition(All2$Label, p = 0.65, list = FALSE, times = 1)
# Training set
train<-Betas_select2[train_ind, ]
# Test set
test<-Betas_select2[-train_ind, ]
xtrain <- train[ , !(names(train) %in% "Label")] %>% as.matrix()
xtest <- test[ , !(names(test) %in% "Label")] %>% as.matrix()
# xtrain<-model.matrix(Label~. -1, data = train)
ytrain<-train$Label
# xtest<-model.matrix(Label~. -1, data = test)
ytest<-test$Label
xtest<-data.frame(xtest)
xtrain<-data.frame(xtrain)
y<-Betas_select2[,NumVar]
X<-Betas_select2[,1:(NumVar-1)]
levels(ytrain)[1]<-"0"
levels(ytrain)[2]<-"1"
levels(ytest)[1]<-"0"
levels(ytest)[2]<-"1"
# source('AUCFun29-3.R')
## multiple, good
multipleAUCNB[1,j]<-multipleAUCfunNB(xtrain, ytrain,xtest,ytest)
multipleAUC[1,j]<-multipleAUCfun(xtrain, ytrain,xtest,ytest)
multipleNBROC[[j]]<-multipleROCfunNB(xtrain, ytrain,xtest,ytest)
multipleROC[[j]]<-multipleROCfun(xtrain, ytrain,xtest,ytest)
#separate
## single, should be good
for (s in (1:(NumVar-1))){
s<-as.numeric(s)
singleAUC[s,j]<-singleAUCfun(xtrain, ytrain,xtest,ytest,s)
singleROC[[s]]<-singleROCfun(xtrain, ytrain,xtest,ytest,s)
MatsinglePlus[s,j]<-singlePlusfun(xtrain, ytrain,xtest,ytest,s)
}
MatsingleROC[,j]<-matrix(singleROC)
#
# for (s in (1:(NumVar-1))){
# cat(s)
# s<-as.numeric(s)
# k<-CompareISS
# doubleAUC[s,j]<-doubleAUCfun(xtrain, ytrain,xtest,ytest,s,k)
# doubleROC[[s]]<-doubleROCfun(xtrain, ytrain,xtest,ytest,s,k)
# MatdoublePlus[s,j]<-doublePlusfun(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCSVM[s,j]<-doubleAUCfunSVM(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCSVMCross[s,j]<-doubleAUCfunSVMCross(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCNB[s,j]<-doubleAUCfunNB(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCRFCross[s,j]<-doubleAUCfunRFCross(xtrain, ytrain,xtest,ytest,s,k)
#
# }
#
#
# MatdoubleROC[,j]<-matrix(doubleROC)
#
#
#Null Hypothesis ####
# Training set
train$Label<-sample(train$Label)
test$Label<-sample(test$Label)
#Permuted data, will make sure that are models are really valid as randomizing the label should yield around 0.5 AUC values. The same testing and training arrangements for the real per model are used.
# Test set
### John Check dim(train) here
xtrain <- train[ , !(names(train) %in% "Label")] %>% as.matrix()
xtest <- test[ , !(names(test) %in% "Label")] %>% as.matrix()
# xtrain<-model.matrix(Label~. -1, data = train)
ytrain<-train$Label
# xtest<-model.matrix(Label~. -1, data = test)
ytest<-test$Label
xtest<-data.frame(xtest)
xtrain<-data.frame(xtrain)
y<-Betas_select2[,NumVar]
X<-Betas_select2[,1:(NumVar-1)]
levels(ytrain)[1]<-"0"
levels(ytrain)[2]<-"1"
levels(ytest)[1]<-"0"
levels(ytest)[2]<-"1"
# source('AUCFun29-3.R')
multipleAUCNBR[1,j]<-multipleAUCfunNB(xtrain, ytrain,xtest,ytest)
multipleAUCR[1,j]<-multipleAUCfun(xtrain, ytrain,xtest,ytest)
multipleNBROCR[[j]]<-multipleROCfunNB(xtrain, ytrain,xtest,ytest)
multipleROCR[[j]]<-multipleROCfun(xtrain, ytrain,xtest,ytest)
for (s in (1:(NumVar-1))){
s<-as.numeric(s)
singleAUCR[s,j]<-singleAUCfun(xtrain, ytrain,xtest,ytest,s)
singleROCR[[s]]<-singleROCfun(xtrain, ytrain,xtest,ytest,s)
MatsinglePlusR[s,j]<-singlePlusfun(xtrain, ytrain,xtest,ytest,s)
}
MatsingleROCR[,j]<-matrix(singleROCR)
#
# for (s in (1:(NumVar-1))){
# s<-as.numeric(s)
# k<-CompareISS
# doubleAUCR[s,j]<-doubleAUCfun(xtrain, ytrain,xtest,ytest,s,k)
# doubleROCR[[s]]<-doubleROCfun(xtrain, ytrain,xtest,ytest,s,k)
# MatdoublePlusR[s,j]<-doublePlusfun(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCSVMR[s,j]<-doubleAUCfunSVM(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCSVMCrossR[s,j]<-doubleAUCfunSVMCross(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCNBR[s,j]<-doubleAUCfunNB(xtrain, ytrain,xtest,ytest,s,k)
# doubleAUCRFCrossR[s,j]<-doubleAUCfunRFCross(xtrain, ytrain,xtest,ytest,s,k)
#
# }
#
#
# MatdoubleROCR[,j]<-matrix(doubleROCR)
#
#
}
##### End Block 5
##### Begin Block 6
#if (NISShere==1){
# names2<-sapply(1:(NumVar-1), function(i){paste0("ISS/",names(xtrain)[i])})
#}else{
names2<-sapply(1:(NumVar-1), function(i){paste0("NISS/",names(xtrain)[i])})
#}
trial<-NULL
h=1
for (b in (2:NumVar-1)) {
print(b)
trial[b]<-names2[h]
h=h+1
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
# doubleAUC<-as.data.frame(doubleAUC)
# doubleAUC<-mutate(doubleAUC, Means=rowMeans(doubleAUC))
# row.names(doubleAUC)<-trial
singleAUC<-as.data.frame(singleAUC)
singleAUC<-mutate(singleAUC, Means=rowMeans(singleAUC))
row.names(singleAUC)<-names(xtrain)
# doubleAUCR<-as.data.frame(doubleAUCR)
# doubleAUCR<-mutate(doubleAUCR, Means=rowMeans(doubleAUCR))
# row.names(doubleAUCR)<-trial
singleAUCR<-as.data.frame(singleAUCR)
singleAUCR<-mutate(singleAUCR, Means=rowMeans(singleAUCR))
row.names(singleAUCR)<-names(xtrain)
multipleAUCNBR<-as.data.frame(multipleAUCNBR)
multipleAUCNBR["Means"]<-rowMeans(multipleAUCNBR)
multipleAUCR<-as.data.frame(multipleAUCR)
multipleAUCR["Means"]<-rowMeans(as.data.frame(multipleAUCR))
multipleAUCNB<-as.data.frame(multipleAUCNB)
multipleAUCNB["Means"]<-rowMeans(as.data.frame(multipleAUCNB))
multipleAUC<-as.data.frame(multipleAUC)
multipleAUC["Means"]<-rowMeans(as.data.frame(multipleAUC))
Final<-data.frame(MultiNB=t(multipleAUCNB),MultiNBRand=t(multipleAUCNBR),Multi=t(multipleAUC),MultiRand=t(multipleAUCR) )
FinalMeans<-data.frame(MultiNB=multipleAUCNB$Means,MultiNBRand=multipleAUCNBR$Means,Multi=multipleAUC$Means,MultiRand=multipleAUCR$Means )
#sacar ROC CURVES #####
print("here")
## [1] "here"
for (g in 1:(NumVar-1)){
plot(MatsingleROC[[g,1]],lwd=3,main=paste("ROC curve of", names(xtrain)[g]))
for (b in 1:N){
plot(MatsingleROCR[[g,b]],col=b,lty=3,add=TRUE)
}
}
#grid.arrange(plot3[[l]], arrangeGrob(plot4[[l]],plot[[l]], ncol=2), ncol=1)
# for (t in 1:(NumVar-1)){
# cat("\n")
# cat(t)
# plot(MatdoubleROC[[t,1]],lwd=3,main=paste("ROC curve of", trial[t]))
# for (b in 1:N){
# plot(MatdoubleROCR[[t,b]],col=b,lty=3,add=TRUE)
# }
#
# }
plot(multipleROC[[1]],lwd=3,main=paste("ROC curve of", trial[1]))
for (b in 1:N){
plot(multipleNBROCR[[b]],col=b,lty=3,add=TRUE)
}
#
# #P value of systematic analysis
#
# max1<-(c(rep(0,(NumVar-1))))
# max2<-(c(rep(0,(NumVar-1)*2)))
#
# for(f in 1:(NumVar-1)){
# for (h in 1:N){
# if (singleAUCR[f,h]>singleAUC$Means[f]){
# max1[f]<-max1[f]+1
# }
# }
# }
#
# for(f in (NumVar-1)){
# for (h in N){
# if (doubleAUCR[f,h]>doubleAUC$Means[f]){
# max2[f]<-max2[f]+1
# }
# }
# }
#
#
#
AUCT<-as.data.frame(t(singleAUC[1:dim(singleAUC)[2]-1]))
AUC2T<-as.data.frame(t(singleAUCR[1:dim(singleAUCR)[2]-1]))
one<-as.data.frame(singleAUC$Means)
two<-as.data.frame(singleAUCR$Means)
Mean<-data.frame(one,two)
Mean<-as.matrix(Mean)
#
nm<-names(AUCT)
for (j in 1:(NumVar-1)){
Mono<-data.frame(Mono=AUCT[,j],Label=as.factor(c(rep("model",dim(AUCT)[1]))))
Mono1<-data.frame(Mono=AUC2T[,j],Label=as.factor(c(rep("permuted data",dim(AUC2T)[1]))))
Mono<-rbind(Mono,Mono1)
print(ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+geom_vline(xintercept =Mean[j,],linetype = 2,color="black",show.legend = TRUE)+labs(title= paste("Permuted versus real model, density plot of AUC curve of",nm[j]),y=paste("Density",nm[j]),x="AUC values")+geom_text(aes(x=as.numeric(Mean[j,1]),y=-0.25),label=signif(Mean[j,1], digits = 2))+geom_text(aes(x=as.numeric(Mean[j,2]),y=-0.25),label=signif(Mean[j,2], digits = 2)))
}
#
# AUCTd<-as.data.frame(t(doubleAUC[1:dim(doubleAUC)[2]-1]))
# AUC2Td<-as.data.frame(t(doubleAUCR[1:dim(doubleAUCR)[2]-1]))
# three<-as.data.frame(doubleAUC$Means)
# four<-as.data.frame(doubleAUCR$Means)
# Meand<-data.frame(three,four)
# Meand<-as.matrix(Meand)
# nmd<-names(AUCTd)
# #
# #
# for (j in 1:(NumVar-1)){
# cat("\n")
# cat(j)
# Monod<-data.frame(Monod=AUCTd[,j],Label=as.factor(c(rep("model",dim(AUCTd)[1]))))
# Mono1d<-data.frame(Monod=AUC2Td[,j],Label=as.factor(c(rep("permuted data",dim(AUC2Td)[1]))))
# Monod<-rbind(Monod,Mono1d)
# print(ggdensity(Monod, x = "Monod", fill = "Label", palette = "jco")+geom_vline(xintercept=Meand[j,],linetype = 2,color="black",show.legend = TRUE)+labs(title= paste("Permuted versus real model, density plot of AUC curve of",nmd[j]),y=paste("Density",nmd[j]),x="AUC values")+geom_text(aes(x=as.numeric(Meand[j,1]),y=-0.25),label=signif(Meand[j,1], digits = 2))+geom_text(aes(x=as.numeric(Meand[j,2]),y=-0.25),label=signif(Meand[j,2], digits = 2)))
# }
# #
n <- length(plot)
# nCol <- floor(sqrt(n))
# # do.call("grid.arrange", grobs=c(plots, ncol=nCol))
# # grid.arrange(grobs = plot, ncol = 2)
# ###Accuracy, Precision etc
#
plot3<-list()
plot4<-list()
Data<-c("Accuracy","Sensitivity","Specificity","Precision")
##### End Block 6
##### Begin Block 7
#
#Accuracy
Mono<-NULL
MonoR<-NULL
for (s in 1:(NumVar-1)){
for (d in 1:4){
for (j in 1:N){
Mono1<-data.frame(Mono=MatsinglePlus[s,j][[1]][[d]])
Mono<-rbind(Mono,Mono1)
Mono1R<-data.frame(MonoR=MatsinglePlusR[s,j][[1]][[d]])
MonoR<-rbind(MonoR,Mono1R)
}
MonoLab<-data.frame(Mono=Mono,Label=as.factor(c(rep("model",dim(Mono)[1]))))
MonoLabR<-data.frame(Mono=MonoR,Label=as.factor(c(rep("permuted data",dim(MonoR)[1]))))
colnames(MonoLabR)[1]<-c("Mono")
Means<-data.frame(Mean=mean(MonoLab$Mono),MeanR=mean(MonoLabR$Mono))
Mono<-rbind(MonoLab,MonoLabR)
plot3[[d]]<-ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+labs(title= paste("Value : ",Data[d], "for", names(xtrain)[s]),x=paste(names(xtrain)[s]))
#+geom_vline(xintercept=Means[1,],linetype = 2,color="black",show.legend = TRUE)+geom_text(aes(x=as.numeric(Means[1,1]),y=-0.25),label=signif(Means[1,1], digits = 2))+geom_text(aes(x=as.numeric(Means[1,2]),y=-0.25,color="blue"),label=signif(Means[1,2], digits = 2))
Mono<-NULL
MonoR<-NULL
}
print(grid.arrange(plot3[[2]], plot3[[3]],plot3[[4]],plot3[[1]], ncol=2))
}
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
#plot3<-ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+labs(title= paste("Value studied single: ",Data[d], "for", names(xtrain)[s]))
#
# Mono<-NULL
# MonoR<-NULL
# for (s in 1:(NumVar-1)){
# for (d in 1:4){
# for (j in 1:N){
# Mono1<-data.frame(Mono=MatdoublePlus[s,j][[1]][[d]])
# Mono<-rbind(Mono,Mono1)
# Mono1R<-data.frame(MonoR=MatdoublePlusR[s,j][[1]][[d]])
# MonoR<-rbind(MonoR,Mono1R)
#
# }
# MonoLab<-data.frame(Mono=Mono,Label=as.factor(c(rep("model",dim(Mono)[1]))))
# MonoLabR<-data.frame(Mono=MonoR,Label=as.factor(c(rep("permuted data",dim(MonoR)[1]))))
# colnames(MonoLabR)[1]<-c("Mono")
# Mono<-rbind(MonoLab,MonoLabR)
# plot4[[d]]<-ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+labs(title= paste("Value : ",Data[d], "for", trial[s]),x=paste(trial[s]))
# Mono<-NULL
# MonoR<-NULL
# }
# print(grid.arrange(plot4[[2]], plot4[[3]],plot4[[4]],plot4[[1]], ncol=2))
#
# }
#
Meanq<-data.frame(multipleAUC$Means,multipleAUCR$Means)
MA<-data.frame(Mono=t(multipleAUC[1:dim(multipleAUC)[2]-1]))
MA["Label"]<-as.factor(c(rep("model",dim(MA)[1])))
MAR<-data.frame(Mono=t(multipleAUCR[1:dim(multipleAUCR)[2]-1]))
MAR["Label"]<-as.factor(c(rep("permuted data",dim(MAR)[1])))
Mono<-rbind(MA,MAR)
pp<-ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+geom_vline(xintercept=Meanq[1,1],linetype = 2,color="black",show.legend = TRUE)+labs(title= paste("AUC curve of GLM Multiple",names1),y="Density",x="AUC values")+geom_text(aes(x=as.numeric(Meanq[1,1]),y=-0.25),label=signif(Meanq[1,1], digits = 2))+geom_text(aes(x=as.numeric(Meanq[1,2]),y=-0.25),label=signif(Meanq[1,2], digits = 2))
print(pp+geom_vline(xintercept=Meanq[1,2],linetype = 2,color="black",show.legend = TRUE))
Meanq<-data.frame(multipleAUCNB$Means,multipleAUCNBR$Means)
MA<-data.frame(Mono=t(multipleAUCNB[1:dim(multipleAUCNB)[2]-1]))
MA["Label"]<-as.factor(c(rep("model",dim(MA)[1])))
MAR<-data.frame(Mono=t(multipleAUCNBR[1:dim(multipleAUCNBR)[2]-1]))
MAR["Label"]<-as.factor(c(rep("permuted data",dim(MAR)[1])))
Mono<-rbind(MA,MAR)
pp<-ggdensity(Mono, x = "Mono", fill = "Label", palette = "jco")+geom_vline(xintercept=Meanq[1,1],linetype = 2,color="black",show.legend = TRUE)+labs(title= paste("AUC curve of NB Multiple",names1),y="Density",x="AUC values")+geom_text(aes(x=as.numeric(Meanq[1,1]),y=-0.25),label=signif(Meanq[1,1], digits = 2))+geom_text(aes(x=as.numeric(Meanq[1,2]),y=-0.25),label=signif(Meanq[1,2], digits = 2))
print(pp+geom_vline(xintercept=Meanq[1,2],linetype = 2,color="black",show.legend = TRUE))
ee<-regex(names,".")
# NamesModels<-c("doubleAUCNB","doubleAUCNBR","doubleAUCRFCross","doubleAUCRFCrossR","doubleAUCSVM","doubleAUCSVMR","doubleAUCSVMCross","doubleAUCSVMCrossR")
# q<-list(doubleAUCNB,doubleAUCNBR,doubleAUCRFCross,doubleAUCRFCrossR,doubleAUCSVM,doubleAUCSVMR,doubleAUCSVMCross,doubleAUCSVMCrossR)
# s<-lapply(q,function(i){MeansNames(i,trial)})
# Models<-data.frame(s[[1]][,dim(doubleAUCNB)[2]])
# for (h in 1:length(s)){
# Models[,h]<-data.frame(s[[h]][,dim(doubleAUCNB)[2]])
# }
# row.names(Models)<-trial
# names(Models)<-NamesModels
Final<-data.frame(MultiNB=t(multipleAUCNB),MultiNBRand=t(multipleAUCNBR),Multi=t(multipleAUC),MultiRand=t(multipleAUCR) )
FinalMeans<-data.frame(MultiNB=multipleAUCNB$Means,MultiNBRand=multipleAUCNBR$Means,Multi=multipleAUC$Means,MultiRand=multipleAUCR$Means )
##### End Block 7
### John code to print Laura's functions
#sink("LauraFunctions.R")
#for (jk in seq_along(lsf.str())) {
# myFun <- lsf.str()[jk]
# cat(paste(myFun," <- ",sep = ""))
# print(getAnywhere(myFun)[1])
#}
#sink()
genes_to_network <- colnames(Betas_select)[-length(colnames(Betas_select))]
genes_to_networkStrict <- colnames(Betas_select2)[-length(colnames(Betas_select2))]
sink("genes_to_network.txt")
cat("Beta_Select_Intial")
## Beta_Select_Intial
cat("\n")
cat(genes_to_network, sep = "\n")
## SULF1
## INHBA
## COL10A1
## COL8A1
## FN1
## NTM
## NOX4
## THBS2
## ADAMTS12
## RASAL2
## CAPG
## FAP
## CTHRC1
## VCAN
## WISP1
## TIMP1
## LTBP1
## SLPI
## GPRC5A
## MIR34AHG
## AEBP1
## KRT7
cat("Beta_Select_Threshold")
## Beta_Select_Threshold
cat("\n")
cat(genes_to_networkStrict, sep = "\n")
## GPRC5A
## CAPG
## INHBA
## COL8A1
## SULF1
## RASAL2
## SLPI
## ADAMTS12
## THBS2
## VCAN
## LTBP1
## AEBP1
## MIR34AHG
## KRT7
sink()